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ABSTRACT 

i 

Theoretically based turbulence models have had success in predicting many features of 
incompressible, free shear layers. However, attempts to extend these models to the high- 
speed, compressible shear layer have been less effective. In the present work, the compressible 
shear layer was studied with a second-order turbulence closure, which initially used only vari- 
able density extensions of incompressible models for the Reynolds stress transport equation 
and the dissipation rate transport equation. The quasi-incompressible closure was unsuc- 
cessful; the predicted effect of the convective Mach number on the shear layer growth rate 
was significantly smaller than that observed in experiments. Having thus confirmed that 
compressibility effects have to be explicitly considered, a new model for the compressible 
dissipation was introduced into the closure. This model is based on a low Mach number, 
asymptotic analysis of the Navier-Stokes equations, and on direct numerical simulations of 
compressible, isotropic turbulence. The use of the new model for the compressible dissipation 
led to good agreement of the computed growth rates with the experimental data. Both the 
computations and the experiments indicate a dramatic reduction in the growth rate when 
the convective Mach number is increased. Experimental data on the normalized maximum 
turbulence intensities and shear stress also show a reduction with increasing Mach number. 
The computed values are in accord with this trend. 


1 This research was supported by the National Aeronautics and Space Administration under NASA Con- 
tract No. NAS1-18605 while the author was in residence at the Institute for Computer Applications in 
Science and Engineering (ICASE), NASA Langley Center, Hampton, VA 23665. 
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1 Introduction 


The reduced growth rate of the high-speed, compressible shear layer relative to its low-speed 
counterpart has been confirmed in several experimental studies, for example, in the recent 
investigations of Papamoschou and Roshko 1 , and Elliott and Samimy 2 . However, variable 
density extensions of incompressible turbulence models, without any explicit compressibility 
terms, have failed to predict the significant decrease in the spreading rate caused by an 
increase in the convective Mach number. This has led to attempts by Oh 3 , Vandromme 4 , 
and Dussauge and Quine 5 , among others, to make phenomenological modifications to in- 
compressible turbulence models, in order to obtain successful predictions of the compressible 
mixing layer. Recently, Sarkar et al. 6 and Zeman 7 have recognized the importance of an addi- 
tional contribution to the turbulent dissipation rate, which is generated by the non-negligible 
fluctuating dilatation in compressible turbulence. The additional term - the compressible dis- 
sipation - has been modeled by Sarkar et al. 6 ; this model is based on a low Mach number, 
asymptotic analysis of the compressible Navier-Stokes equations and is calibrated with ref- 
erence to direct numerical simulations of compressible, isotropic turbulence. The present 
paper applies the model of the compressible dissipation to the high-speed shear layer within 
the framework of a second-order turbulence closure. A schematic of the shear layer is given 
in Fig. 1. 

The paper is organized in the following manner. In Section 2 the exact governing equa- 
tions are given, and the turbulence models constituting the second-order closure are de- 
scribed. The numerical procedure is outlined in Section 3. The results of the calculations 
with the second-order closure are given in Section 4, and conclusions are presented in Sec- 
tion 5. 

2 The governing equations 

We obtain the equations for the mean variables by first decomposing each variable into a 
mean component and a fluctuating component, and then averaging the equations for the 
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following variables: the density p, the velocity Ui and the total energy E. The total energy 
E is defined by 

E = -^ + C v T (1) 

where T denotes the static temperature, and C v is the specific heat at constant volume. 
The Reynolds decomposition of an instantaneous variable f> into its mean and fluctuating 
components is 

4> = t + <t>" 

where, by definition, <f>" = 0. The Favre decomposition of an instantaneous variable is also 
used in compressible turbulence, primarily because the resulting structure of the averaged 
inertial terms is simpler; this decomposition is given by 


4> = l + <t>' 

where <f> is the density-weighted Reynolds average, 

P 

The overbar over a variable is used to denote a conventional Reynolds average, while the 
overtilde is used to denote the Favre average. A single superscript ' represents fluctuations 
with respect to the Favre average, while a double superscript " signifies fluctuations with 
respect to the Reynolds average. The conventional Reynolds average of Favre fluctuations 


is non-zero, in particular, <f>' = —p"4>"lp. After averaging the instantaneous Navier-Stokes 
equations, the following mean equations are obtained: 

Conservation of mass: 

dt(p) + (puk),k = 0 (2) 

Conservation of momentum: 


d t (pui ) + (pujfcUi),* = -p ti + T ikik - (pu-O ifc 
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where the anisotropy tensor bij is given by 


_ ufoj Sjj 

XJ q 2 3 

and q 2 = denotes the trace of the Reynolds stress tensor. In (9) the model 

constants are 

Ci — 3.0 , = 0.6 

Since the primary aim of the paper is to study the influence of terms that arise solely from 
flow compressibility, we do not use more sophisticated incompressible pressure-strain models, 
such as those proposed by Shih and Lumley 11 ; Fu, Launder and Tselepidakis 12 ; and Speziale, 
Sarkar and Gatski 13 . 

The dissipation rate tensor tij is commonly believed to be isotropic at high turbulence 
Reynolds numbers, leading to the model 

€ij = ( 10 ) 
where the turbulent dissipation rate e is given by 


pe — T kl u kl 


( 11 ) 


The viscous stress in a compressible fluid is 

2 

T ij = — '^f J ' U k,k^ij (1^) 

where we have neglected the bulk viscosity. As shown in Sarkar et al. 6 , substitution of (12) 
into (11), followed by some algebraic manipulation, gives 

pe = p(e, + e c ) (13) 


where 


e, — vw"u}" 


and 

e c — \vd" 2 
3 


(14) 

(15) 
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Here a is the fluctuating vorticity, and d" = u kk is the fluctuating dilatation. The decom- 
position (13) of the turbulent dissipation rate e into the solenoidal dissipation e 4 and the 
compressible dissipation e c is asymptotically valid for high- Reynolds number turbulence, and 
is exact for constant- viscosity homogeneous turbulence. Because of the explicit compressible 
contribution to the turbulent dissipation rate, the treatment of e has to be modified with 
respect to the incompressible case. Developing an appropriate, direct modification of the 
transport equation for e is a difficult proposition, because the exact transport equation for e 
is complicated for the incompressible case, and even more so for the compressible case. Also, 
as discussed by Speziale 14 , the addition of new terms into the e transport equation has often 
led to unintended, deleterious effects in homogeneous flow. In the present work, we adopt 
a simpler alternative. The incompressible form of the dissipation equation is retained as a 
transport equation for e„; such an approach is valid, because e, is not affected by moderate 
levels of compressibility 6 . It remains to model e c ; we choose the simple, algebraic model of 
Sarkar et al. 6 , 

e c = a x e,M f ( 16 ) 

which is motivated by an asymptotic analysis of the compressible Navier-Stokes equations 
with M t as the small parameter. Here M t denotes the turbulent Mach number defined by 
M t = \Jq 2 /-yRT, and f is the Favre-averaged temperature. Finally, the model for e i5 becomes 

2 

e ij = 3 P£* ( 1 + M] )6ij (17) 

The model constant was set as orj = 1 with reference to direct numerical simulations of the 
decay of isotropic, compressible turbulence. Zeman 7 has also used a similar decomposition 
of the turbulent dissipation rate, and after assuming that eddy shocklets occur in high-speed 
flows, he derives a model for the contribution of these eddy shocklets to the compressible 
dissipation. 

In the present work, we assume that theFuIk viscosity p v = 0. If the bulk viscosity p v is 
non-negligible, for example in polyatomic gases, there is an additional turbulent dissipation 
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term = JT^d" 2 which can be modeled as ej, = a 2 e a M 2 . If the value of //„ is known, a 2 can 
be easily determined from a x by the relation a 2 = 3/vafi/4/i. 

The pressure-dilatation p"d", which is not necessarily single-signed (i.e; it is neither 
positive semi-definite or negative semi-definite) like the compressible dissipation, is a more 


difficult term to model. Low Mach number asymptotic theory 15,6 suggests that p"d" is 
negligible compared to e c , and from direct simulations 6 it appears that in isotropic, moderate 
Mach number turbulence p"d" is appreciably smaller than e c . In the present closure, we will 
neglect p"d" relative to e c . 

The diffusive transport T^k is modeled by a gradient transport expression, 

M 2 ) 2 - - - 


T ijk = -c,p — [0'u') ifc + («x),i + OX)J 


where C, — 0.018. The quantity u[ is related to the turbulent mass flux p"u" by 


u' = 


P"< 


and after using (6) for the mass flux, we obtain the model 


c„k 2 _ 

- — p.i 

peG p 


I ^ U 

u ~ 


(18) 


(19) 


( 20 ) 


The standard high-Reynolds number form of the dissipation rate equation is used as the 
transport equation for e 5 , 

dt{pz.) + {pUke,),k = -C € i ^-putyu^ - C t2 p y + (C e — u'^e, ) >k (21) 

The model coefficients in (21) are 

C €l = 1.44 , C e2 = 1.90 , C e = 0.15 (22) 


For the present problem, we need to solve (2)-(4), along with the equation of state, 
to obtain the mean variables: ~p, U , V , and E. In the case of the plane shear layer, the 
Reynolds stress tensor has four non-zero components: u'v\ u' 2 , v' 2 and w' 2 , which are solved 
by the corresponding components of (8). The equation for the solenoidal dissipation rate e 3 
completes the set of governing equations. Thus a system of nine coupled, non-linear, partial 
differential equations along with an appropriate set of initial and boundary conditions must 
be solved. 
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3 Method of Solution of the Governing Equations 


The transport equations for the mean flow and Reynolds stresses are written in the phys- 
ical domain and must be transformed to the computational domain using an appropriate 
coordinate transformation. For the physical problem under consideration, an algebraic grid 
generation technique is used to generate the mesh. In the physical domain a uniform grid 
is used in the axial direction and in the normal direction the grid lines are clustered near 
regions where strong gradients exist. A uniform mesh is used in the computational domain. 
The governing equations are first cast into a vector form, where U is the dependent variable 
vector consisting of nine components, the vectors F and G respectively denote the x and y 
flux vectors, and H is the source vector containing the terms causing production, destruc- 
tion and redistribution of the Reynolds stresses. To numerically obtain the solution for the 
vector i7, the governing equations are then transformed from the physical domain to the 
computational domain, giving the following system of equations, 


dU dF dG « 

aT + W + W 


(23) 


where 


U = JU , H = JH 

F = Fy v - Gx v , G = Gx i - Fy ( , J = x£y n - y^x v . 

In (23), a superscript Q denotes quantities in the transformed system, (x^, x nJ y^ } y r7 ) represent 
the metrics of the transformation, and J denotes the Jacobian of the transformation. If the 
physical grid is given, the metrics and the Jacobian of the transformation can be easily 
computed. 

The governing equations are integrated explicitly in time using the unsplit MacCormack 
predictor- corrector scheme. During a specific numerical sweep, the inviscid fluxes and the 
first-derivative terms in the source vector H are backward differenced in the predictor step 
and forward differenced in the corrector step. Second-order central differences are used for 
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the viscous and heat flux terms. Hence the complete scheme for both the predictor and 

corrector steps can be expressed as follows 

Predictor. 


ATT " +1 _ Aj f | V T]GiJ T 

AU " - - At {-Kr + -zr~ u 

T ? n+1 TT 71 , A TT n+1 

Uij — Uij + AUij 


Corrector. 


A K 


n+l 


= -At 


AA n+i , a 


+ 


A£ ' A77 

* n+l 1 / ; n * n+1 A * n+l\ 

Uij = ~ J 


- 


n+1 


The composite numerical scheme is second-order accurate in both time and space and, 
being an explicit scheme, is conditionally restricted by the Courant and viscous stability 
limits of the governing equations. The solution procedure requires no scalar or block tridi- 
agonal inversions. The flow field is advanced from time level n to n+1 and this process is 
continued until the desired integration time or steady state has been reached. Since the 
Reynolds stress transport equations contain stiff source terms, the maximum CFL number 
used in the computation was limited to 0.5. 

The numerical code used in this study is a two-dimensional, elliptic, Navier-Stokes solver 
(SPARK2D 16 ) written in a generalized body-oriented coordinate system. As such, various 
two-dimensional free shear flows and wall bounded flows can be handled by the numerical 
code. The code in its original form used a second-order spatially and temporally accurate, 
two-step MacCormack scheme. The latter versions of the code employ a variety of higher- 
order compact algorithms 17 (4th and 6th order) and various upwind schemes. Local time 
stepping and residual smoothing options are also available in the code to accelerate the 
convergence to steady state. Both laminar reacting and non-reacting flows can be easily 
handled by the code. In the present research work, the capabilities of SPARK2D are further 
enhanced by adding a second-order Reynolds stress model as a turbulence closure. 


9 



Since the governing equations are elliptic in nature, the boundary conditions have to be 
specified along all four boundaries. These include inflow, outflow and outer boundaries (lower 
and upper boundaries) respectively. At the inflow boundary (x=0.0), profiles are specified 
for the velocities, static pressure, static temperature, turbulent stresses and the turbulent 
dissipation rate. Since we are interested in the downstream fully-developed regime, the 
specific form of the inlet profiles is not crucial. 

The outer boundaries always remain in the free-stream and the appropriate boundary 
condition is to assume that the normal derivative of the flow variables vanish along those 
boundaries. The gradient boundary conditions, not only preserve the free-stream values 
along the outer boundaries but also provide nonreflective conditions for the outgoing waves. 
The outflow boundary ( x = x max ) is always supersonic, and hence the values of mean flow and 
turbulence quantities are determined by zeroth-order extrapolation from upstream values. 
Along with the boundary conditions, the governing equations also require a set of initial 
conditions. The initial conditions at time t=0 for all the variables are obtained by simply 
propagating the inflow profiles throughout the computational domain. Having specified all 
the boundary and initial data the equations are marched in time until the residual based 
on ~pU decreases by six orders of magnitude, indicating that a converged solution has been 
obtained. 


4 Results 


It is known that the fully-developed, high- Reynolds number shear layer spreads linearly, and 
that the growth rate dS/dx satisfies the relation 


d ± = C ( U '~ U ^ 
dx n U 1 + U 2 ) 


(24) 


where 6(x ) denotes the width of the shear layer, and C £ is approximately constant. The 
shear layer thickness S(x ) has been defined in several ways by previous investigators; in the 
present work, £(z) represents the distance between the two cross-stream positions where the 
normalized streamwise velocity U* = (U — U 2 )/(Ui — U 2 ) is respectively 0.1 and 0.9. The 
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fully-developed nature of the shear layer is also characterized by the maximum values of the 
normalized turbulent stresses <7 U , <j v , cr w and cr uv reaching constants; where 



c 7 UV = \j —U 9 V 9 !(U\ - U 2 ) 


Figs. 2-6 show results for a particular set of conditions for the shear layer between two 
streams of air. The high-speed stream had a velocity lh = 2500 m/s while the low-speed 
stream had a velocity U 2 = 800 m/s. The thermodynamic quantities in the two incident 
streams were equal and were prescribed as T\ = 800 K, pi = 1 atm, and pi = 0.44 kg/m . 
When the ratio of specific heats 7 has the same value in the two streams, the convective 
Mach number M c is given by 1 , 


Mr = 


Ui-U 2 


0,1 + a 2 

where a\ and a 2 are the respective speeds of sound in the two layers. The case described by 
Figs. 2-6 corresponds to M c — 1.5. The computational domain for this case was a rectangle 
of dimensions 0.1 m x 0.05 m with a 201 x 51 grid overlaying it. The grid spacing was uniform 
in the streamwise direction and stretched in the cross-stream direction. Based on comparison 
with results using other grid spacings, the resolution of the 201 X 51 grid for the computational 
domain was found to be sufficient to provide practically grid-independent results for the mean 
velocity and turbulent stress profiles. As an example of the grid sensitivity of the calculated 
solution, increasing the number of grid points by a factor of approximately 1.7 changed the 
values of C&, and the maximum values of cr u , cr v , a w and cr uv by less than 2 % from the values 
corresponding to the 201 x 51 grid. 

Fig. 2 shows that the shear layer thickness S(x ) increases linearly after an initial develop- 
ment phase. In Fig. 3 the normalized streamwise mean velocity U* at the inlet, outlet and 
two intermediate locations is plotted as a function of the similarity variable 77 = (y — y c )/£> 
where y is the local cross-stream coordinate and y c is the cross-stream coordinate where 
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U* = 0.5. It Is evident from Figs. 2 and 3 that, at the outflow boundary of the computa- 
tional box, the linearly growing regime is well-established and the mean velocity has reached 
its self-similar profile. The similarity mean velocity profile of Fig. 3 is somewhat asymmetric 
with respect to its center rj = 0 and indicates a greater penetration into the low-speed side 
than the corresponding penetration into the upper, high-speed side of the domain. Fig. 4 
shows the mean temperature profile across the shear layer. There is a sharp increase of the 
temperature in the core of the shear layer due to the large velocity gradients there. Figs. 5 
and 6 show profiles of the normalized streamwise turbulence intensity a u and the normalized 
shear stress cr uv . All the components of the normalized Reynolds stress tensor reach their 
self-similar profiles at the exit of the computational box. 

The growth rate parameter C& and the maximum values of the normalized Reynolds 
stresses a uy cr v , a w and a uv are nominally constant for the incompressible shear layer. How- 
ever, it is clear from the experimental data of Figs. 7 and 8 that these quantities show a 
systematic decrease when the convective Mach number M c increases. In Fig. 7, the incom- 
pressible value (Cj) o, which was obtained by calculating a case with a small M c , was used 
to normalize the growth rate parameter C$. Fig. 7 indicates that the Reynolds stress calcu- 
lations without the compressibility model (16) show only a modest decrease in the growth 
rate parameter. However, introduction of the model for the compressible dissipation leads 
to good agreement with both the experimentally observed trends of the sharp decrease in 
the growth rate, and the later flattening of the growth rate curve in the high Mach num- 
ber range. It is evident from Fig. 8 that computations with the compressible dissipation 
model are in qualitative agreement with the observed trend of a decrease in the maximum 
normalized Reynolds stress components with an increase in M c . 

Growth rate curves for various values of are shown in conjunction with the Langley 
experimental curve 18 in Fig. 9. Increasing ati from its recommended value of 1.0 leads to a 
sharper reduction of the growth rate before the eventual flattening out at high convective 
Mach numbers. The flattening of the growth rate curve for high M c is due to the maximum 
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turbulent Mach number M t asymptoting to an equilibrium level (as shown in Fig. 10), and 
consequent leveling out of the compressible contribution to the turbulent dissipation rate. 

The model of Sarkar et al. 6 for the compressible dissipation, which was used in the 
present work, has also been applied by Wilcox 19 to some supersonic and hypersonic flows 
within the framework of a k — u turbulence closure. Wilcox’s study concludes that the 
addition of this model of the compressible dissipation leads to the experimentally observed 
reduction in the growth rate of the compressible shear layer, leads to values of skin friction in 
adiabatic boundary layers that are somewhat lower than the measured values, and results in 
an improved prediction of the separation bubble size in a shock-boundary layer interaction 
problem. 

5 Conclusions 

Initially, a second-order turbulence closure without any explicit compressibility models was 
applied to the high-speed shear layer. The results confirmed earlier conclusions 18,20,21 re- 
garding the inability of such variable density generalizations of incompressible models to 
predict the strong influence of the convective Mach number on the growth rate of the shear 
layer. The new model of Sarkar et al. 6 for the compressible dissipation was then incorpo- 
rated into a full Reynolds stress closure. The growth rates computed with this model, not 
only captured the experimentally observed sharp reduction of the growth rate at interme- 
diate Mach numbers, but also showed the tendency to flatten out at large Mach numbers. 
The present calculations are also in agreement with the experimental result that the maxi- 
mum normalized turbulence intensities and shear stress decrease when the convective Mach 
number is increased. 

In the future, we propose to apply the present second-order closure to more complex 
compressible flows. Though, the consequences of the enhanced dissipation in compressible 
flows are consistent with some of the distinguishing features of the high-speed shear layer, 
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other compressibility phenomena may become important in different flows like the shock- 
boundary layer interaction. Our future studies will address issues relevant to the modeling 
of such distinct mechanisms. 
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Figure 2. Downstream evolution of the shear layer thickness. 
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Figure 4. Transverse mean temperature profiles at various streamwise locations 




Figure 5. Transverse profiles of the streamwise component of the Reynolds stress tensor. 
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Figure 6. Transverse profiles of the Reynolds shear stress. 
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Figure 7. Variation of the growth rate of the compressible shear layer with the convective 
Mach number. 
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Figure 9. Computed growth rate curves for various values of the parameter ai in the 
model for compressible dissipation. 
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